function stability_analysis()
    % 设置步长数组
    h_values = [0.2, 0.1, 0.05, 0.001];
    % 问题区间
    x_range = [0, 1];
    % 初值
    y0 = 1;
    
    % 创建图形窗口
    figure('Position', [100, 100, 1000, 600]);
    
    % 计算精确解用于比较
    x_exact = linspace(0, 1, 1000);
    y_exact = exp(-20 * x_exact);
    
    % 对每个步长计算并绘图
    for i = 1:length(h_values)
        h = h_values(i);
        
        % 构建x网格
        x = x_range(1):h:x_range(2);
        
        % Euler方法
        y_euler = euler_method(x, y0, h);
        
        % 4阶龙格库塔方法
        y_rk4 = rk4_method(x, y0, h);
        
        % 创建子图
        subplot(2, 2, i);
        plot(x, y_euler, 'b.-', 'DisplayName', 'Euler方法');
        hold on;
        plot(x, y_rk4, 'r.-', 'DisplayName', 'RK4方法');
        plot(x_exact, y_exact, 'k-', 'DisplayName', '精确解');
        
        title(sprintf('h = %.3f', h));
        xlabel('x');
        ylabel('y');
        grid on;
        legend('Location', 'best');
        
        % 计算并显示误差
        euler_error = max(abs(y_euler - exp(-20*x)));
        rk4_error = max(abs(y_rk4 - exp(-20*x)));
        fprintf('\n步长 h = %.3f:\n', h);
        fprintf('Euler方法最大误差: %.2e\n', euler_error);
        fprintf('RK4方法最大误差: %.2e\n', rk4_error);
    end
end

function y = euler_method(x, y0, h)
    % Euler方法
    n = length(x);
    y = zeros(1, n);
    y(1) = y0;
    
    for i = 1:(n-1)
        y(i+1) = y(i) + h * (-20*y(i));
    end
end

function y = rk4_method(x, y0, h)
    % 4阶龙格库塔方法
    n = length(x);
    y = zeros(1, n);
    y(1) = y0;
    
    for i = 1:(n-1)
        k1 = -20 * y(i);
        k2 = -20 * (y(i) + h*k1/2);
        k3 = -20 * (y(i) + h*k2/2);
        k4 = -20 * (y(i) + h*k3);
        
        y(i+1) = y(i) + h * (k1 + 2*k2 + 2*k3 + k4) / 6;
    end
end